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Detailed computations of tokamak edge turbulence in three dimensional, 
globally consistent flux tube geometry show an inhibition of the standard 
scenario in which zonal ExB flows generated by the turbulence should lead to 
transport barrier formation. It is found by comparison to slab geometry and 
by analysis of the energetics that the zonal flow energy is depleted by toroidal 
coupling to the pressure through the geodesic curvature. Edge transport 
barriers would then depend on the physics of the neoclassical equilibrium. 
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Drift Wave Turbulence and Zonal Flows. Drift wave turbulence is nonlinear, nonperiodic 
motion involving disturbances on a background thermal gradient of a magnetised plasma 
and eddies of fluid like motion in which the advecting velocity of all charged species is the 
ExB velocity [1,2]. The disturbances in the electric field implied by the presence of these 
eddies are caused by the tendency of the electron dynamics to establish a force balance 
along the magnetic field. Pressure disturbances have their parallel gradients balanced by a 
parallel electric field, whose static part is given by the parallel gradient of the electrostatic 
potential. This potential in turn is the stream function for the ExB velocity in drift 
planes, which are locally perpendicular to the magnetic field. The turbulence is driven by 
the background gradient, and the electron pressure and electrostatic potential are coupled 
together through parallel currents. Departures from the static force balance are mediated 
primarily through electromagnetic induction and resistive friction, but also the electron 
inertia, which is not negligible [3] . Further details are provided by the temperatures, whose 
dynamics is very robust due to nonlinear, time dependent Landau damping. In a three 
dimensional, toroidal flux surface geometry, the turbulence is characterised by a nonlinear 
instability whose inherent vorticity is strong enough to "supersede" linear interchange 
instabilities, giving tokamak edge turbulence a drift wave basic character [4]. 

Although this turbulence effects an unsteady transport through nonlinear advection of 
the thermodynamic state variables, actual modification of the profiles proceeds on the much 
slower transport time scale, typically at least three orders of magnitude slower than the 
turbulence even in steep gradient regions. Such quasilinear modification of the background 
(a three wave interaction involving a wave and its complex conjugate driving changes in 
the background) also occurs in the other variables, specifically the ExB vorticity, for which 
the time scale is short enough for a self-consistent interaction with the turbulence to affect 
the dynamics of the turbulence itself. These are the "zonal flows" [5] , which are simply the 
ExB flows resulting from disturbances in the electrostatic potential which are constant on 
a given magnetic flux surface. One can think of a zonal flow as a rigid poloidal rotation 
of the entire flux surface. These flows are important because when and where they are 
sheared they can cause a local suppression of the turbulence. Suppression of turbulence by 
sheared flows began by considering an imposed flow which is part of the background [6] . It 
was then pointed out that the process by which this suppression occurs conserves energy, 
mainly involving a transfer of energy in three wave interactions from smaller scale eddies 
and the larger scale background [7], a variant of the more general inverse energy cascade 
from smaller to larger scales in two dimensional, incompressible turbulence [8] . Suppression 
by imposed ExB shear was then shown in computations to proceed energetically [9]. Self 
consistent interactions concern time dependent zonal flows which have time and space 
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scales comparable to or only somewhat larger than those of the turbulence. Their study as 
such [10], followed global scale computational studies showing them to be very important 
in limiting the radial scale of the turbulence and consequently the resulting transport [11]. 
Their importance in fusion research lies in the fact that zonal ExB flow (electric field) 
shear is thought to underly the transition and maintenance of the H-mode operation of 
tokamak confinement [12]. 

We here examine the physics of the zonal flow/ turbulence interaction in a model of the 
turbulence which includes two important generalisations: departures from the "adiabatic" 
state of perfect electron force balance, and an electromagnetic character in that response, 
which allows significant delays in the adiabatic response at larger perpendicular scales, 
since the collisional response is through the parallel current, while the inductive response 
is through the time dependence of the parallel magnetic potential (by Ampere's law, the 
current is given by the perpendicular Laplacian of that potential). We will find that in 
toroidal geometry, although the drive of the zonal flows via Reynolds stress remains, the 
geodesic curvature of the magnetic field lines couples the zonal flows to pressure sidebands 
with finite parallel gradient. These sidebands serve as a localised part of the general source 
for the turbulence, as the free energy transfer in the pressure disturbances is preferentially 
towards smaller scales. The build up of strong, long lived ExB "mean flow" shear layers 
is thereby inhibited, preventing the turbulence from self consistently generating enough 
ExB shear to strongly reduce its own amplitude. This prevents the scenario in which the 
zonal flow drive process should lead to transport barrier formation. In slab geometry, the 
geodesic curvature effect is absent and mean flows do develop, but we find by inserting 
specifically this geodesic curvature effect that the toroidal result is recovered (incidentally 
demonstrating the weakness of the ballooning/interchange effect in the turbulence). It is 
important to note in this context that models of transport barrier formation by Reynolds 
stress-induced self-generated flows rely on two-dimensional slab geometry [13], and they 
work well in such computations [14], but the three-dimensional toroidal result is rather 
different as documented herein. 

The DALF3 Model. The simplest three dimensional model of drift wave dynamics which 
takes the self consistent adiabatic response into account and allows it an electromagnetic 
character is a four field model in toroidal flux tube geometry called DALF3 (the drift Alfven 
model [3] but omitting the temperature dynamics). The state variables are the electrostatic 
potential and the electron pressure p e , and the flux variables are the parallel current Jy 
and the parallel ion velocity u\\ , all expressed as disturbances on the equilibrium which is 
a set of constant parameters except where the ExB and magnetic nonlinearities operate on 
the background gradients. In the Ohm's law, electromagnetic induction, electron inertia, 



3 



and resistive friction are all retained. The adiabatic response is the reaction of the parallel 
current, controlled by those three effects, to the pressure/potential static force imbalance, 
acting to couple p e and through the shear Alfven dynamics. The model equations in 
simplified sheared flux tube geometry are 

n e MiC 2 d „ 9 ~ J\\ „ c „ . . 

1 dAy m e dJu ~ 1 _ ~ 

I (Pe = y£V,i -PeBV„^ - iv • ^BxVp e +p e V • ^BxV? (3) 

n e Mi-^- = -V,| (p e +p e ) (4) 

with Ampere's law Jy = — (c/47r)Vj_A||. The ExB advective and parallel derivatives are 
given by 

S = S + ^ BXV *' V V^iB.V-^Bxviln-V (5) 

where B and B are the equilibrium magnetic field and its magnitude, and the combinations 
involving and Ay give the ExB and disturbed parallel derivatives, respectively, i.e., the 
nonlinear ities. The flux tube geometry used is detailed elsewhere [15], as is the importance 
of global consistency which controls field line connection [16]. The standard normalisation 
is in terms of the drift scale p s = c s /Vti and frequency c s /Lj_, where c 2 s = T e /Mi and 
fij = eB/MiC, and L± is the background scale length for p e . The parameters controlling 
the adiabatic response are f3 = (c s /L±) 2 (qR/vA) 2 , and jl = (c s /L±) 2 (qR/V e ) 2 , and C = 
0.51(z/ e Lj_/c s )/i, reflecting the competition between perpendicular ExB turbulence and the 
parallel dynamics, where the field line connection length is 2nqR, V e is the electron thermal 
velocity (V 2 = T e /m e ) and the 0.51 comes from the parallel resistivity, r]\\ = 0.51m e z/ e /n e e 2 
[17]. The sound waves are controlled by e = (c s / L±) 2 (qR/c s ) 2 , just the parallel/perp scale 
ratio. The effects of magnetic curvature (the radius of curvature is R, the toroidal major 
radius), entering through /C = V- (c/B 2 )BxV are controlled by ujb = 2L±/R, which 
can be set independently - - slab geometry is ujb = 0. The coordinates are {x, y, s}, 
representing the down- gradient, electron drift, and parallel directions, respectively. The 
computations are set up exactly as detailed in [15], with a grid of 64 x 256 x 16 nodes 
in {x, y, s}, and with node spacings h x = h y = 207rp s /64 and h s = 2nqR/lQ. Nominal 
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parameters corresponding to a typical plasma edge in the L-mode of tokamak operation 
are 

(3 = 2 fl = 5 C = 7.65 u B = 0.05 e = 18350 s = 1 (6) 
roughly reflecting physical parameters: 

n e = 4.5 x 10 13 cm -3 T e = 80 eV S = 2.5T (7) 

i? = 165 cm L± = 3.65 cm q = 3 (8) 

This is rather strongly collisional (standard parameter = 40), but because Cub < 1 it 
is still well within the drift wave regime [18]. 
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Figure 1. Time traces of the squared amplitudes of <fi (Ap), p e (An), and Vj_</> 
(Aw). Due to the disparate k± factors, Ap tracks mostly the flows and Aw 
mostly the turbulence. The basic slab case shows initial saturation and then 
weakening of the turbulence as the flow amplitude rises. The basic toroidal 
case shows persistent saturation, as the flow amplitude remains low. 

We refer to the cases with ujb = and 0.05 as the basic slab and toroidal cases, 
respectively. The time traces for these cases appear in Fig. 1. The squared amplitudes 
are shown for (Ap), p e (An), and the vorticity Vj_0 (Aw). When strong zonal flow 
layers appear, they dominate the Ap signal because of the lack of k± factors. The Aw 
signal by contrast, with four additional k± factors, mostly tracks the turbulence. For the 
basic slab case, the turbulence saturates in the time range 200 < t < 400, after which it 
is ground down by the rise of the overall flow levels; Ap grows to large values, and Aw 
correspondingly decreases. For the basic toroidal case, the saturation occurs at roughly 
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the same time scale, but the Ap curve saturates unsteadily at a much lower level, smaller 
by about two orders of magnitude as in the slab case, reflecting the flows which are simply 
part of the turbulence. The turbulence saturates and maintains its level, close to the basic 
gyro Bohm transport. All time traces reflect this saturated state. We therefore find that 
the spin up and suppress scenario operates moderately well in slab geometry but not at 
all in toroidal geometry (for the same basic result in models including both temperatures 
see [19]). 

The effort to explain this perhaps startling result forces systematic address of the 
various toroidal effects, all of which (in this model) operate through the curvature terms. 
There are two basic effects in the curvature operator /C: the interchange dynamics itself, 
and the geodesic curvature. Pure interchange dynamics operates on the k y 7^ part, 
through IC y d/dy. The geodesic curvature effect is in IC x d/dx. To test directly for the 
geodesic curvature one must separate K, x out for the k y = part and leave the IC y pure 
interchange effect alone. The reason for suspecting the geodesic curvature is that the zonal 
flow effects lie in the k y = part, for which the pure interchange effect vanishes due to the 
vanishing d/dy. 

The Geodesic Curvature Effect. The basic mode of oscillation involving the geodesic 
curvature is the classic MHD geodesic acoustic oscillation [20], which represents simple 
coupling between the pressure and vorticity through the geodesic curvature K. x . The 
pressure part of this is a sideband with parallel wavenumber k\\qR = 1 which presents 
itself as an interim free energy source to the turbulence. The potential part is the zonal 
flow; both are axisymmetric (k y = 0). Due to the strong direct cascade tendency, the 
nonlinear ExB pressure advection, we ■ Vp e , quickly delivers this free energy back to the 
turbulence. Overall, this transfer process acts as a depletion channel for zonal flow energy, 
keeping the zonal flow amplitude at levels comparable to the turblence. The loss channel 
is from the zonal flow to the sideband p e (the IC X terms in Eqs. 1,3), and then through 
ve • Vp e back to the eddies of the turbulence. The mutual energy transfer is conservative, 
so the tendency of the system to reach equipartition results in a finite population of the 
zonal flow mode, but not so large as to overwhelm the turbulence. We note that the fact 
that the geodesic curvature couples all fey modes of a given k y requires us to keep or remove 
it for all k y = modes as a unit; otherwise, the resulting model would not conserve energy. 

A modified toroidal case is constructed by taking IC X out of the k y = part of the 
basic toroidal case, and a modified slab case is made by putting K, x into the k y = part 
of the basic slab case, thereby isolating the geodesic curvature effects on the k y = part. 
The results concerning the turbulence amplitudes and transport are shown in Fig. 2. Two 
time traces are shown for each case: the transport (Qe) and the <fi squared amplitude (Ap). 
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Figure 2. Time traces of the squared amplitudes of <fi (Ap) and the transport 
(Qe), for the four cases. The basic toroidal and modified slab cases reach 
persistent saturation; both contain the geodesic curvature effect. The basic 
slab and modified toroidal cases lack this effect and are both ground down by 
strong, self generated flow shear. This test confirms the geodesic curvature 
effect of coupling zonal flows to finite ku sidebands as the reason the spin up 
and suppress scenario does not work in toroidal geometry. 

We find immediately that the two cases without geodesic curvature in the k y = part are 
similar, with the flow amplitude rising to high values, grinding down the transport. The 
two cases with the geodesic curvature in the k y = part are also similar, with the potential 
amplitude kept at levels low enough that the turbulence is not suppressed. 

The morphology of <p and p e is shown for the four cases in Fig. 3. The basic slab 
case shows dominance of the shear layer in the potential, with the vorticity disturbances 
stretched into thin sheets sharply tilted into the y-direction. The basic toroidal case 
also shows shear layers, but their vorticity represents a frequency not larger than that of 
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Figure 3. Morphology of the flows and disturbances (1/2 of the y-domain is 
shown). The basic slab case shows the strong shear layers in 0, and sheets 
of vorticity (Vj_</>) stretched in the y-direction. The basic toroidal case shows 
visible shear layers in (p, but they are of similar magnitude as the turbulence 
and do not strongly affect the form of V^_</>. These are the time dependent 
zonal flow layers visible in the unsteady <p amplitude in Fig. 2. The modified 
toroidal case appears slablike, while the modified slab case looks like the basic 
toroidal case. 
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the basic turbulence, which is why a strong amount of suppression does not occur. The 
modified toroidal case shows the strong shear layer of the basic slab case, and the modified 
slab case shows the structure of the basic toroidal case. The shear levels of these weaker 
flows are comparable to the dynamical frequencies of the turbulence (about 0.1c s /L±). 
The Ap curves for these two cases with weaker flows show those flows to be short-lived, 
comparable to the correlation time of the turbulence (about 6L±/c s ). These are the zonal 
flows which remain as part of the turbulence, leading in fact to moderate suppression but 
allowing it to remain at a robust amplitude. 
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Figure 4. Snapshots of the zonal flow profiles for the four cases. The basic 
toroidal and modified slab cases show the weak, time dependent zonal flows 
which are part of the turbulence. The basic slab and modified toroidal cases 
show the strong, self generated shear layers which suppress the turbulence. 



The instantaneous profiles (the zonal flow mode, k y = fey = 0) of 4> are shown for the 
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four cases in Fig. 4, in the same arrangement as for the time traces. The flow shear of 
these is strong or weak according to whether the geodesic curvature is absent or present 
in the k y = part, respectively. A rough guide of whether these sheared flow layers are 
able to suppress the turbulence is given by what can be called the "diamagnetic flow shear 
level" given by 

n D = v D /L ± (9) 

where vp = cT e /eBLj_ is the diamagnetic velocity. The level of shear in the 4> profiles 
(actually given by the vorticity profile) is well below this for the two unsuppressed cases 
(basic toroidal and modified slab), and well above this for the other two cases. (ALT: show 
vor profiles, note diag shear is p* in n.u., while turb omega is about 0.1) 
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Figure 5. Transport trend versus collisionality. The drift wave regime extends 
to v = 10, at which Cojb ~ 1, and for these parameters the standard is 
136. The basic toroidal and slab cases are compared to companion runs in 
which the zonal flow drive is removed. The time dependent zonal flows are the 
difference in the toroidal case; the self generated shear layers, in the slab case. 
The comparison between the basic toroidal and modified slab cases shows the 
role of pure interchange forcing; the difference at v = 3 (C = 7.65) is about 20 
percent. 



The transport results for various collisionality are summarised in Fig. 5 (note C = 
2.55z/). The basic toroidal and slab cases show similar trends if zonal flows are eliminated 
entirely by removing the flux surface average of we • VV^</>. The fluctuating zonal flows 
provide a slightly reduced transport in the basic toroidal case, but the strong shear layers 
in the basic slab case strongly suppress the turbulence, even more so for larger C. The 
modified slab case is much like the basic toroidal case, showing the effects of geodesic 
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curvature to inhibit the strong shear layers, leaving the fluctuating zonal flows and the 
pure interchange effects intact. The small difference between those two cases is the pure 
interchange effect, incidentally showing that direct interchange drive in toroidal geometry 
is but a small perturbation on an existing drift wave mode structure. 

The coupling of zonal flows to geodesic acoustic oscillations can be further demon- 
strated through the energy theorem satisfied by the zonal flows and the pressure sidebands. 
Let (• • •) denote the flux surface average and note that it commutes with d/dx but an- 
nihilates d/dy. Let {• • •} further denote the average over x. The zonal flow potential is 
(0), the zonal flow is (v y ) = (d<fi/dx), the zonal vorticity is (O) = {dv y /dx) : and the zonal 
flow energy is j(i> y ) 2 j. Through the geodesic curvature the zonal flow is coupled to the 
Pfirsch-Schliiter Alfven mode and then again to modifications in the background pressure 
(assuming unit diagonal metric, and neglecting sound waves, magnetic nonlinear ities, and 
sidebands with |&||<Z-R| > 1): 

Wt {I {vV)2 } = m ^ vV ^- ub {^ sins > ( 10 ) 

{(p e sins) 2 } = 2 j^^sins^ (<2 x sins) J + u B {(p e sins) (v y )} 

-u B |(p e sins) ^-t^ j ~ 2 {(Pe sins)(J|| coss)} 



(11) 



d_ 
di 



d_ 
di 



j/T 1 (B y coss) 2 + fi(J\\ coss) 2 } 

= 2 {(J\\ cos s) ((p e — 4>) sins)} — 2C |( J\\ cos s) 2 | 

! } = 2 {(O sins) (v x v y sms)} 

+ ujb |^^^^ (0 sins) I + 2 {(</>sins) (Jy coss)} 



(12) 



(v y sin s) 



(13) 



8 r WW^W-^%>W-Pe)sin S H (14) 



dt {2 J { \ dx / J [ \ dx 

where Q x = p e v x is the pointwise transport and B y — —f3dA\\/dx is the field disturbance. 
If time averages are taken, the left sides of these equations vanish and the right sides become 
balances between drive, transfer, and depletion mechanisms. The drive for the zonal flow 
is the zonally averaged Reynolds stress (v x v y ), correlated with the zonal vorticity. The 
depletion mechanism is geodesic transfer to the Pfirsch-Schliiter sideband (p e sins), and 
the depletion for that is the nonlinear transfer of free energy back to the turbulence (which 
requires nonadiabatic electrons, enabled by the finite (3). The flow sideband (v y sins) is 



11 



controlled by relaxation of the Pfirsch-Schluter current (Jycoss). The proper geodesic 
acoustic oscillation is that between ((f)) and {p e sin s) as correctly noted in [20] , the only 
subsystem not strongly affected by the Alfven dynamics. We recognise {(dp e /dx) (Q x )} 
as the negative of the main drive of the self sustained turbulence (cf. [21]), and hence 
as the corresponding depletion of the background (p e ). The profile is maintained by the 
damping of the k y = part of p e to zero at the boundaries in x, which affects both 
(p e ) and (p e sms). For the nominal case the zonal Reynolds stress and geodesic transfer 
were measured at 0.858 ± 0.359 and 1.02 ± 0.305, and depletion of the sideband went 
through the nonlinearity, the profile maintenance, and the Pfirsch-Schluter transfer to 
the background at 0.470 ± 0.431 and 0.301 ± 0.100 and 0.103 ± 0.0741, respectively, with 
all other sideband effects much smaller (all numbers xlO -2 ). The nonlinear depletion 
of (p e sins) functions because of the pointwise correlation of Q x with —dp e /dx. Indeed, 
while either of these two sin s terms is small in time average, the average of their product is 
not, because of this correlation. This energetic depletion therefore overwhelms any slight 
presence of a nonzero (<5 x sins) (Stringer- Taylor effect). Indeed, the PDF of (<5 x sins) 
was found to be close to Gaussian. The pressure nonlinearity is therefore a depletion 
of the sideband and therefore ultimately of the zonal flow. The energy flow is p e — > 
(<fi) — > (p e sins) — > p e , through the Reynolds stress, geodesic curvature, and nonlinear flux 
correlation, respectively. The Stringer- Taylor effect is therefore a sink for the zonal flow 
system, not a source as incorrectly reported in Ref. [22], whose runs were apparently not 
taken to complete statistical saturation and in any case suffer from all the shortcomings 
of the drift resistive ballooning paradigm (cf. Ref. [18]). 

Main Points. The principal result of this study is that while the turbulent Reynolds 
stress always tends towards a transfer of energy from small eddies to the larger scale 
zonal flows (similar k x but disparate k y ), in toroidal geometry the geodesic curvature 
couples the zonal flows to finite-fey pressure sidebands, which act as a loss channel by 
means of nonlinear advective transfer back to the turbulence. This prevents large scale, 
large amplitude zonal flows from forming and therefore rules out the spin up and suppress 
scenario for transport barrier formation, at least due to local (homogeneous) action of the 
ExB Reynolds stress. The ExB shear layers observed in tokamak edge transport barriers 
[12] must therefore come from some other mechanism, most likely having to do with the 
neoclassical equilibrium. Two recent proposals are a generalised ion orbit loss mechanism 
[23,24], and the generation of a large parallel flow and its equilibrium electric field profile 
by coupling to the open field line regions [25]. For the core regions the electron response 
is more adiabatic and electrostatic, so that current results on core zonal flows [11] are not 
affected. 
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